clear
cap log close
set more off

/***********************************************************************************
***	This program examines rotation group bias in the Current Population Survey   ***
*** using samples linked over time as presented in                               ***
*** Krueger, Mas, and Niu (Review of Economics and Statistics, forthcoming).     ***
************************************************************************************/

// Change to the directory that contains datasets
local datadir "C:"
// Change to the directory that contains output files
local outdir "C:"

// Log file
log using `outdir'/RGB_CPS_codes_linked_sample_REStat.log, replace

// Set see for bootstrapped SE
set seed 1234567890

/***********************************************************************************
	Merge individuals over time
************************************************************************************/

// Append annual survey files across years

#delimit ;
local vars94 year intmonth monis hhid1 hhid2 hhid04 hhnum pulineno
	         age female race unemp lf weight
			 intstthh noninta94;
local vars76 year intmonth monis hhid1 hhid2 hhnum pulineno
	         age female race unemp lf weight;
#delimit cr

foreach y of numlist 1994/2014 {
	di "Year `y'"
	use `vars94' using `datadir'/base`y'ax.dta, clear
	gen intstat = 1 if (pulineno < .)	// Residing in an interviewed households
	compress
	save `outdir'/temp`y'.dta, replace
}
foreach y of numlist 1976/1993 {
	di "Year `y'"
	use `vars76' using `datadir'/base`y'ax.dta, clear
	gen intstat = 1 if (pulineno < .)
	compress
	save `outdir'/temp`y'.dta, replace
}

// Append annual data
use `outdir'/temp1976.dta, clear
foreach y of numlist 1977/2014 {
	append using `outdir'/temp`y'.dta	
}
gen month = ym(year,intmonth)
replace hhid04 = 0 if year <= 1993
save `outdir'/temp.dta, replace

// Erase annual data
foreach y of numlist 1976/2014 {
	erase `outdir'/temp`y'.dta
}

// Extract data from each rotation group (MIS 1-8)
foreach i of numlist 1/8 {

	* Keep one rotation group (include those not interviewed)
	use `outdir'/temp.dta if monis == `i', clear
	drop monis
	
	* Calculate the month first selected into the survey	
	if `i' <= 4 {
		replace month = month-`i'+1
	}
	else {
		replace month = month-`i'+1-8
	}
	
	* Year and month of entry into CPS
	gen starty = year(dofm(month))
	gen startm = month(dofm(month))
	drop month
	
	// Rename variables other than household and person identifiers (attaching MIS number)
	rename age      age`i'
	rename female   female`i'
	rename race     race`i'
	rename intstat  intstat`i'	
	rename unemp unemp`i'
	rename lf    lf`i'	
	rename weight weight`i'
	
	rename intstthh intstthh`i'
	rename noninta94 noninta94`i'
	
	#delimit ;
	keep starty startm hhid1 hhid2 hhid04 hhnum pulineno
		age`i' female`i' race`i' intstat`i' unemp`i' lf`i' weight`i'
		intstthh`i' noninta94`i';
	#delimit cr
	
	* Sort data by ID
	sort starty startm hhid1 hhid2 hhid04 hhnum pulineno
	
	compress
	save `outdir'/temp`i'.dta, replace
	
	// The sample of interviewed individuals
	use `outdir'/temp`i'.dta, clear
	keep if intstat`i' == 1
	* Remove duplicate observations
	quietly by starty startm hhid1 hhid2 hhid04 hhnum pulineno: gen dup = cond(_N==1,0,_n)
	tab dup
	tab starty dup if dup != 0
	keep if dup == 0	
	drop dup	
	save `outdir'/temp1`i'.dta, replace
	
	// The sample of non-interviewed households
	use `outdir'/temp`i'.dta, clear
	keep if intstat`i' != 1
	* Remove duplicate observations
	quietly by starty startm hhid1 hhid2 hhid04 hhnum: gen dup = cond(_N==1,0,_n)
	tab dup
	tab starty dup if dup != 0
	keep if dup == 0
	drop dup
	save `outdir'/temp2`i'.dta, replace
	
	erase `outdir'/temp`i'.dta
	
}

erase `outdir'/temp.dta

* Merge data across rotation groups
* Interviewed (1:1 merge)
use `outdir'/temp11.dta, clear
merge 1:1 starty startm hhid1 hhid2 hhid04 hhnum pulineno using `outdir'/temp12.dta, nogen
merge 1:1 starty startm hhid1 hhid2 hhid04 hhnum pulineno using `outdir'/temp13.dta, nogen
merge 1:1 starty startm hhid1 hhid2 hhid04 hhnum pulineno using `outdir'/temp14.dta, nogen
merge 1:1 starty startm hhid1 hhid2 hhid04 hhnum pulineno using `outdir'/temp15.dta, nogen
merge 1:1 starty startm hhid1 hhid2 hhid04 hhnum pulineno using `outdir'/temp16.dta, nogen
merge 1:1 starty startm hhid1 hhid2 hhid04 hhnum pulineno using `outdir'/temp17.dta, nogen
merge 1:1 starty startm hhid1 hhid2 hhid04 hhnum pulineno using `outdir'/temp18.dta, nogen

* Not interviewed (m:1 merge)
merge m:1 starty startm hhid1 hhid2 hhid04 hhnum using `outdir'/temp21.dta, nogen update
merge m:1 starty startm hhid1 hhid2 hhid04 hhnum using `outdir'/temp22.dta, nogen update
merge m:1 starty startm hhid1 hhid2 hhid04 hhnum using `outdir'/temp23.dta, nogen update
merge m:1 starty startm hhid1 hhid2 hhid04 hhnum using `outdir'/temp24.dta, nogen update
merge m:1 starty startm hhid1 hhid2 hhid04 hhnum using `outdir'/temp25.dta, nogen update
merge m:1 starty startm hhid1 hhid2 hhid04 hhnum using `outdir'/temp26.dta, nogen update
merge m:1 starty startm hhid1 hhid2 hhid04 hhnum using `outdir'/temp27.dta, nogen update
merge m:1 starty startm hhid1 hhid2 hhid04 hhnum using `outdir'/temp28.dta, nogen update

foreach i of numlist 1/8 {
	erase `outdir'/temp1`i'.dta
	erase `outdir'/temp2`i'.dta
}

// Indicator variables for presence in surveys
gen in1 = (intstat1 == 1)
gen in2 = (intstat2 == 1)
gen in3 = (intstat3 == 1)
gen in4 = (intstat4 == 1)
gen in5 = (intstat5 == 1)
gen in6 = (intstat6 == 1)
gen in7 = (intstat7 == 1)
gen in8 = (intstat8 == 1)
* Present in all eight interviews
gen in1_8 = (in1 & in2 & in3 & in4 & in5 & in6 & in7 & in8)

// Find observations with consistently reported demographic infomration (Age and gender)
egen maxage1_8 = rowmax(age1 age2 age3 age4 age5 age6 age7 age8)
egen minage1_8 = rowmin(age1 age2 age3 age4 age5 age6 age7 age8)

egen maxfemale1_8 = rowmax(female1 female2 female3 female4 female5 female6 female7 female8)
egen minfemale1_8 = rowmin(female1 female2 female3 female4 female5 female6 female7 female8)

gen age1_8    = (maxage1_8 - minage1_8   <= 2 | maxage1_8 == .)
gen female1_8 = (maxfemale1_8 == minfemale1_8 | maxfemale1_8 == .)

// Age and gender are consistent
gen SA1_8 = (age1_8 & female1_8)

// In all interviews and demographic information matches
gen in1_8_demo = (in1_8 & SA1_8)

* Drop unused variables
drop min*1_8 max*1_8 age1_8 female1_8

compress
save `outdir'/RGB_CPS_link_interviews.dta, replace

/***********************************************************************************
	Figure A1: annual match rate
************************************************************************************/

* set begining and ending data month
local y0 = 1976
local m0 = 1
local y1 = 2014
local m1 = 12

foreach i of numlist 1/8 {
	use `outdir'/RGB_CPS_link_interviews.dta, clear
	
	* Keep those who started the interview from Jan 1976 and 16 months before the last interview month
	keep if ym(starty,startm) >= ym(`y0',`m0') & ym(starty,startm) <= ym(`y1',`m1')-16+1
	* Calculate match rate
	collapse (sum) in`i' in1_8 [pw=weight`i'], by(starty startm)
	gen group = `i'
	gen rate1_8 = in1_8/in`i'

	keep starty startm group rate1_8	
	compress
	save `outdir'/a`i'.dta, replace
}

// Append survey samples across MIS
use `outdir'/a1.dta, clear
foreach i of numlist 2/8 {
	append using `outdir'/a`i'.dta
}

sort starty startm group
order starty startm group rate1_8, first
compress
save `outdir'/RGB_CPS_link_interviews_match_rate_m.dta, replace

// Erase monthly data
foreach i of numlist 1/8 {
	erase `outdir'/a`i'.dta
}

// Plot match rates
* Average match rate by year and MIS
collapse rate1_8, by(starty group)
rename rate1_8 rate
reshape wide rate, i(starty) j(group)
* Average match rate over MIS by year
egen rate_avg = rowmean(rate1-rate8)

// Plot annual match rates
sort starty
#delimit ;
scatter rate1 rate8 rate_avg starty, connect(l l l) sort clpattern(dash dash solid)
	xtitle("Year in which households are first selected into the sample")
	ytitle("Fraction of respondents who are present in all 8 interviews",size(small))
	xlabel(1975(5)2015)
	legend(order(1 "MIS 1" 2 "MIS 8" 3 "Average MIS 1-8") rows(1))
	saving(`outdir'/RGB_CPS_linked_match_rate.gph, replace);
#delimit cr
graph export `outdir'/RGB_CPS_linked_match_rate.pdf, replace

/***********************************************************************************
	Remove survey months in which the match rate is low due to sampling issues
************************************************************************************/
// Find starting month in which the average matching rate is too low (sampling issue rather than match issues)
use `outdir'/RGB_CPS_link_interviews_match_rate_m.dta, clear
* Averate match rate over MIS
collapse rate1_8, by(starty startm)

// If the average match rate across all 8 group is less than 40%, then consider as low
gen low_match1_8 = (rate1_8 < 0.4)

// remove 1994 sample
replace low_match1_8 = 1 if starty == 1994

/* remove pre-1994 cohorts that have observations in 1994+ period
   set begining and ending data month (starting in 1982) */
local y0 = 1982
local m0 = 1
local y1 = 2014
local m1 = 12
#delimit ;
replace low_match1_8 = 1 
	if !((ym(starty,startm) >= ym(`y0',`m0') & ym(starty,startm) <= ym(1993,12)-16+1) |
         (ym(starty,startm) >= ym(1994,1) & ym(starty,startm) <= ym(`y1',`m1')-16+1));
#delimit cr

save `outdir'/low_match1_8.dta, replace

// List sample in the data
list starty startm if low_match1_8 == 0

use `outdir'/RGB_CPS_link_interviews.dta, clear
merge m:1 starty startm using `outdir'/low_match1_8.dta, nogen
drop if low_match1_8 == 1 | low_match1_8 == .
save `outdir'/RGB_CPS_link_interviews_matched.dta, replace

// Delete large data file
erase `outdir'/RGB_CPS_link_interviews.dta

/***********************************************************************************
	Table 4 and Figure B4: Bias and non-interview by demographic group (1994+)
************************************************************************************/

use `outdir'/RGB_CPS_link_interviews_matched.dta if starty >= 1994, clear

foreach i of numlist 1/8 {
	replace race`i' = 3 if race`i' > 2 & race`i' < .
}

// Impute demographic information and sampling weights for non-interviews
foreach i of numlist 1/8 {
	foreach j of numlist 1/8 {
		replace age`i' = age`j'       if weight`i' == . & (weight`j' > 0 & weight`j' < .) & SA1_8
		replace female`i' = female`j' if weight`i' == . & (weight`j' > 0 & weight`j' < .) & SA1_8
		replace race`i' = race`j'     if weight`i' == . & (weight`j' > 0 & weight`j' < .) & SA1_8
		* Update sampling weights
		replace weight`i' = weight`j' if weight`i' == . & (weight`j' > 0 & weight`j' < .) & SA1_8
	}
}
* Get a consistent set of demographic information (the first available non-missing value)
gen age = .
gen female = .
gen race = .
foreach i of numlist 1/8 {
	replace age    = age`i'    if age    == . & age`i'    < .
	replace female = female`i' if female == . & female`i' < .
	replace race   = race`i'   if race   == . & race`i'   < .
}

// Get interview status
foreach i of numlist 1/8 {
	gen interview`i' = (in`i' == 1)
	gen hhtypeA`i' = (intstthh`i' == 2 & !interview`i')
	gen refusal`i' = (noninta94`i' == 3 & hhtypeA`i')	
}

// Keep observations with complete demographic information
recode age (16/24=1) (25/34=2) (35/44=3) (45/54=4) (55/64=5) (65/max=6), gen(agecat)
keep if agecat >= 1 & agecat <= 5	// drop those 65 and older
keep if race == 1 | race == 2
keep if female == 0 | female == 1

save `outdir'/RGB_CPS_linked_demo.dta, replace

**** Non-response rate ****

// Calculate total number of type-A non-interviewed person
foreach i of numlist 1/8 {

	use `outdir'/RGB_CPS_linked_demo.dta, clear
	
	collapse (sum) interview`i' hhtypeA`i' refusal`i' if SA1_8 [pw=weight`i'], by(agecat female race)
	
	// Calculate rate of non-interview
	gen typeArate`i' = hhtypeA`i'/(interview`i' + hhtypeA`i')*100
	gen refusalrate`i' = refusal`i'/(interview`i' + hhtypeA`i')*100
	gen Aother`i' = typeArate`i'-refusalrate`i'
	
	keep agecat female race typeArate`i' refusalrate`i' Aother`i'
	
	compress	
	save `outdir'/temp`i'.dta, replace
	
}

// Merge data across MIS
use `outdir'/temp1.dta, replace
foreach i of numlist 2/8 {
	merge 1:1 agecat race female using `outdir'/temp`i'.dta, nogen
}
* Erase data by MIS
foreach i of numlist 1/8 {
	erase `outdir'/temp`i'.dta
}

* Calculate average rate of non-interview across MIS
reshape long typeArate refusalrate Aother, i(agecat female race) j(monis)
collapse typeArate refusalrate Aother, by(agecat female race)

save `outdir'/temp_nonresp.dta, replace

**** Bias ****

// Calculate the magnitude of RGB
foreach i of numlist 1/8 {
	use `outdir'/RGB_CPS_linked_demo.dta, clear
	collapse (sum) unemp`i' lf`i' [pw=weight`i'] if in`i' & SA1_8, by(agecat female race)
	gen unemprt`i' = unemp`i'/lf`i'
	keep agecat female race unemprt`i'
	save `outdir'/temp`i'.dta, replace
}
// Merge data across MIS
use `outdir'/temp1.dta, replace
foreach i of numlist 2/8 {
	merge 1:1 agecat female race using `outdir'/temp`i'.dta, nogen
}
* Erase data by MIS
foreach i of numlist 1/8 {
	erase `outdir'/temp`i'.dta
}

// Calculate RGB indices of unemployment rate
egen avg = rowmean(unemprt1-unemprt8)
foreach i of numlist 1/8 {
	gen indices`i' = unemprt`i'/avg*100
}

keep agecat female race indices*

// Calculate magnitude of RGB for unemployment rate
reshape long indices, i(agecat female race) j(monis)
gen slope = .
foreach i of numlist 1/5 {
	foreach j of numlist 0/1 {
		foreach k of numlist 1/2 {		
			count if agecat == `i' & female == `j' & race == `k'
			if (r(N)>0) {
				quietly reg indices monis if agecat == `i' & female == `j' & race == `k'
				replace slope = _b[monis] if agecat == `i' & female == `j' & race == `k'
			}
		}
	}
}

collapse slope, by(agecat female race)
save `outdir'/temp_bias.dta, replace

// Merge bias and non-response rates
use `outdir'/temp_nonresp.dta, clear
merge 1:1 agecat female race using `outdir'/temp_bias.dta, nogen

compress
save `outdir'/RGB_CPS_indices_linked_demo.dta, replace

erase `outdir'/temp_nonresp.dta
erase `outdir'/temp_bias.dta

// Table 4: Regression analysis
reg slope typeArate
outreg2 using RGB_CPS_bias_typeA_demo, excel replace
reg slope refusalrate
outreg2 using RGB_CPS_bias_typeA_demo, excel
reg slope Aother
outreg2 using RGB_CPS_bias_typeA_demo, excel

// Figure B4: bias and non-interview rate
#delimit ;
twoway (scatter slope typeArate) (lfit slope typeArate),
	legend(off)
	xtitle("Type-A non-interview rate")
	ytitle("Magnitude of bias")
	title("(a) Type A")
	saving(`outdir'/bias_typeArate, replace);
#delimit cr

#delimit ;
twoway (scatter slope refusalrate) (lfit slope refusalrate),
	legend(off)
	xtitle("Refusal rate")
	ytitle("Magnitude of bias")
	title("(b) Refusal")
	saving(`outdir'/bias_refusal, replace);
#delimit cr

#delimit ;
twoway (scatter slope Aother) (lfit slope Aother),
	legend(off)
	xtitle("Type-A other than refusal")
	ytitle("Magnitude of bias")
	title("(c) Other")
	saving(`outdir'/bias_Aother, replace);
#delimit cr

#delimit ;
graph combine `outdir'/bias_typeArate.gph
              `outdir'/bias_refusal.gph
	          `outdir'/bias_Aother.gph,
			  cols(1) iscale(0.7) imargin(0 0 0 0)
			  xsize(8) ysize(13.5)
			  saving(`outdir'/RGB_CPS_bias_typeA_demo.gph, replace);
#delimit cr
graph export `outdir'/RGB_CPS_bias_typeA_demo.pdf, replace

erase `outdir'/bias_typeArate.gph
erase `outdir'/bias_refusal.gph
erase `outdir'/bias_Aother.gph


/***********************************************************************************
	Table 5: Multiplicative indices by interview participation patterns
************************************************************************************/
* Calculate the magnitude of RGB for 3 groups
* (1) in for all eight interviews
* (2) Other
* (3) All
foreach i of numlist 1/8 {

	* By interview participation pattern
	use `outdir'/RGB_CPS_link_interviews_matched.dta, clear
	collapse (sum) unemp`i' lf`i' [pw=weight`i'] if in`i' & SA1_8, by(starty in1_8_demo)
	gen group = .
	replace group = 1 if in1_8_demo == 1
	replace group = 2 if in1_8_demo == 0
	save `outdir'/temp.dta, replace
	
	* All
	use `outdir'/RGB_CPS_link_interviews_matched.dta, clear
	collapse (sum) unemp`i' lf`i' [pw=weight`i'] if in`i' & SA1_8, by(starty)  // becnchmark everyone
	gen group = 3
	
	* Combine the all groups
	append using `outdir'/temp.dta
	sort starty group
	gen unemprt`i' = unemp`i'/lf`i'
	keep starty group unemprt`i'
	save `outdir'/temp`i'.dta, replace
	
	erase `outdir'/temp.dta	
}

// Merge data across MIS
use `outdir'/temp1.dta
foreach i of numlist 2/8 {
	merge 1:1 starty group using `outdir'/temp`i'.dta, nogen
}
* Erase data by MIS
foreach i of numlist 1/8 {
	erase `outdir'/temp`i'.dta
}

// Calculate indices of RGB for unemployment rate
egen avg = rowmean(unemprt1-unemprt8)
foreach i of numlist 1/8 {
	gen indices`i' = unemprt`i'/avg*100
}
keep starty group indices* avg

// Estimating slope with respect to MIS
reshape long indices, i(starty group avg) j(monis)
gen slope = .
quietly sum starty
local minyr = r(min)
local maxyr = r(max)
foreach y of numlist `minyr'/`maxyr' {
	foreach i of numlist 1/3 {
		count if starty == `y' & group == `i'
		if (r(N)>0) {
			quietly reg indices monis if starty == `y' & group == `i'
			replace slope = _b[monis] if starty == `y' & group == `i'
		}
	}
}

reshape wide indices, i(starty group avg slope) j(monis)

// Average over pre and post-1994 period
gen period = 1 if starty < 1994
replace period = 2 if starty >= 1994
collapse indices1-indices8 avg slope, by(period group)

// Organize output files
order period group indices1-indices8 avg slope

label define period 1 "1982-1993" 2 "1994-2014"
label values period period

label define group 1 "Present MIS1-8" 2 "Missing at least one" 3 "All"
label value group group

compress
save `outdir'/RGB_CPS_indices_linked.dta, replace

log close

/***********************************************************************************
	Table 5: Bootstrapped SE
************************************************************************************/

// Define parameters
local iteration = 500
use `outdir'/RGB_CPS_link_interviews_matched.dta, clear
quietly sum starty
local minyr = r(min)
local maxyr = r(max)

// Start the dataset to store bias by response pattern
clear
set obs 1
gen starty = .
gen iteration = .
gen group = .
gen slope = .
save `outdir'/RGB_CPS_indices_linked_bootstrap.dta, replace
local obscnt = 1

foreach i of numlist 1/`iteration' {

	foreach y of numlist `minyr'/`maxyr' {

		// Sample without weighting (weights change over MIS)
		use `outdir'/RGB_CPS_link_interviews_matched.dta if starty == `y', clear
		
		// Some years have no observatoin betcause of linkage issues
		count
		if (r(N) == 0) {
			continue
		}
		
		// keep if in at least one interview
		keep if in1 | in2 | in3 | in4 | in5 | in6 | in7 | in8
		
		gen wt = .
		bsample, weight(wt)
		
		// Update sampling weights
		foreach j of numlist 1/8 {
			gen wt`j' = weight`j'*wt
			drop weight`j'
		}
		
		compress
		save `outdir'/temp_linked.dta, replace
		
		if (`y' >= 1994) {
			save `outdir'/temp`y'.dta, replace
		}		
		
		foreach j of numlist 1/8 {

			* By interview participation pattern
			use `outdir'/temp_linked.dta, clear
			collapse (sum) unemp`j' lf`j' [pw=wt`j'] if in`j' & SA1_8, by(in1_8_demo)
			gen group = .
			replace group = 1 if in1_8_demo == 1
			replace group = 2 if in1_8_demo == 0
			save `outdir'/temp_local.dta, replace
	
			* All
			use `outdir'/temp_linked.dta, clear
			collapse (sum) unemp`j' lf`j' [pw=wt`j'] if in`j' & SA1_8  // becnchmark everyone
			gen group = 3
	
			* Combine the all groups
			append using `outdir'/temp_local.dta
			sort group
			gen unemprt`j' = unemp`j'/lf`j'
			keep group unemprt`j'
			save `outdir'/temp`j'.dta, replace
	
			erase `outdir'/temp_local.dta	
		}

		// Merge data across MIS
		use `outdir'/temp1.dta
		foreach j of numlist 2/8 {
			merge 1:1 group using `outdir'/temp`j'.dta, nogen
		}
		* Erase data by MIS
		foreach j of numlist 1/8 {
			erase `outdir'/temp`j'.dta
		}

		// Calculate indices of RGB for unemployment rate
		egen avg = rowmean(unemprt1-unemprt8)
		foreach j of numlist 1/8 {
			gen indices`j' = unemprt`j'/avg*100
		}
		keep group indices* avg

		// Estimating slope with respect to MIS
		reshape long indices, i(group avg) j(monis)
		gen slope = .
		foreach j of numlist 1/3 {
			count if group == `j'
			if (r(N)>0) {
				quietly reg indices monis if group == `j'
				replace slope = _b[monis] if group == `j'
			}
		}
		
		collapse slope, by(group)
		save `outdir'/temp_slope.dta, replace
		
		// Record bootstrapped values
		foreach j of numlist 1/3 {
		
			use `outdir'/temp_slope.dta if group == `j', clear
			quietly sum slope
			
			use `outdir'/RGB_CPS_indices_linked_bootstrap.dta, clear
			replace starty = `y' in `obscnt'
			replace iteration = `i' in `obscnt'
			replace group = `j' in `obscnt'
			replace slope = r(mean) in `obscnt'
			
			// Increment the number observations
			local obscnt = `obscnt' + 1
			set obs `obscnt'
			save `outdir'/RGB_CPS_indices_linked_bootstrap.dta, replace
			
		}	
		
		erase `outdir'/temp_slope.dta
		erase `outdir'/temp_linked.dta
	
	}
	
}

log using `outdir'/RGB_CPS_codes_linked_sample_REStat.log, append

// Output files

use `outdir'/RGB_CPS_indices_linked_bootstrap.dta, clear
outsheet using `outdir'/RGB_CPS_indices_linked_bootstrap.csv, comma replace

gen period = 1 if starty < 1994
replace period = 2 if starty >= 1994
keep if starty < .
collapse slope, by(period group iteration)

collapse slope (sd) slope_se = slope, by(period group)
rename slope slope_mean
merge 1:1 period group using `outdir'/RGB_CPS_indices_linked.dta, nogen
order slope_mean slope_se, last
save `outdir'/RGB_CPS_indices_linked.dta, replace
outsheet using `outdir'/RGB_CPS_indices_linked.csv, comma replace


log close
set more on
